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We investigate glassy dynamical properties of one component three-dimensional system of particles 
interacting via pair repulsive potential by the molecular dynamic simulation in the wide region of 
densities. The glass state is superfragile and it has high glassforming ability. The glass transition 
temperature Tg has pronounced minimum at densities where the frustration is maximal. 



The ubiquitous glass formation and jamming still puz- 
zle physicists [3 [5] . The microscopic mechanism of the 
drastic slowing down of the structural relaxation of a liq- 
uid on cooling is one of the central issues of the physics 
of the liquid-glass transition. The question "why some 
liquids form a glass easily but others do not," is still the 
matter of debates. 

There is a paradigm that one component liquid with 
isotropic potential typically spontaneously crystallizes 
being supercooled in (quasi) equilibrium conditions "3'- 
[B]. It is a formidable challenge to avoid, e.g., sponta- 
neous crystallization in quasiequilibrium cooling of the 
onecomponent Lennard- Jones liquid. Yet it was discov- 
ered not long ago that there are some exceptions from 
the paradigm. The common fitch of these exceptions is 
the pronounced attractive well of the pair potential, see, 
e.g., [3|8]. Here we show using the molecular dynamics 
(MD) that the one component simple liquid with pure 
repulsive potential shows glassy behaviour in quasiequi- 
librium cooling. 

Frustration, when one cannot minimize the energy of 
the system by merely minimizing all local interactions, is 
one of the basic factors that stipulates the glass-forming 
ability [Hl^. E.g., it can be related with the long-range 
alternating interactions (e.g., in spin glasses [lOj ) or with 
geometrical reasons [11] . The potential used in Ref. [H |8] 
was optimized to produce icosahedral local order and so 
geometrical frustration. Our potential has the soft step 
and the simple liquid with this interaction has two char- 
acteristic scales. The competition between these scales 
makes the system effectively quasibinary jT5]. And this 
is the origin of frustration in our system. Intuitively in- 
creasing frustration one should favor the formation of the 
glass. Here we show that on the contrary, the glass tran- 
sition temperature Tg may have minimum at parameters 
where the frustration is maximal. 

The second important concept of glassy physics is 
fragility. According to Angell-classification [2 [T3] the 
glassforming liquids effectively divide into two classes: 
"strong" and "fragile", where the viscosity of the liquid 
shows either nearly Arrhenius behaviour with tempera- 
ture or much faster one. We found the fragility index and 
concluded that our system is superfragile [its fragility in- 
dex exceeds that of the "Decalin" - one the most fragile 



liquid [g], see Fig.[T] 

We use the pair potential model of "collapsing soft 
spheres" [ia[ISllin|: 

U{r)=ei^-j +enF[2ko{r-<Ji)], (1) 

where nF{x) = 1/[1 -|-exp(a;)], e - is the unit of energy, a 
and CTi are "hard" -core and "soft" -core diameters. This 
kind of potentials, Eq.Q, is succesfully used for simu- 
lation of water-like anomalies, liquid-liquid phase transi- 
tions and glass formation [121 fT5l - l20] . The graph of the 
potential we discuss in the Supplementary and in Fig. [5] 
In the remainder of this paper we use the dimension- 
less quantities: f = r/cr, U = U/e, temperature T — T/e, 
density p = Ncr^/V, and time i = tl\a ^frnf^, where m 
and V are the molecular mass and system volume corre- 
spondingly. As we will only use these reduced variables, 
we omit the tildes. So we take here n = 14, fcg = 10, and 
CTi = 1.35. These parameters values reveal complex sys- 
tem behavior such as phase diagram with polymorphous 
transitions and disordered gap, see Fig. (2] and water-like 
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FIG. 1. (Color online) Relative fragility of our system as 
the function of density. For comparison we show the fragility 
of the typical glassformers. Our system in the glass regime 
appears to be extremely fragile. It seems that it is in the 
short list of the most fragile glassy systems. 
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FIG. 2. (Color online) We show the glass transition tempera- 
ture on the sketch of the phase diagram obtained in Ref. [T^] • 
The red circles and blue triangles correspond to Tg extracted 
from D{T) using the Vogel-Fulcher formula (VF) [25] and 
the power law from mode-coupling theory (MCT) [31] cor- 
respondingly. In the inset we show the accuracy of D{T) 
approximation by VF and MCT for p — 0.6. The circles show 
the result of MD simulations. 



anomalies P^ITCI ITB]. 

For MD simulations we have used the system of = 
5000 particles that were simulated under periodic bound- 
ary conditions in Nose-Hoover NVT ensemble. We have 
checked that N = 5000 is enough amount of particles to 
eliminate the finite size effects that agrees with Ref. dU 
The MD time step was St — 0.01. It is nearly the maxi- 
mum possible time step that satisfies the energy conser- 
vation condition. The system was studied in the density 
region of p € (0.35 — 0.75). At all densities of this re- 
gion the system was cooled in a stepwise manner from 
high temperature state and completely equilibrated at 
each step until convergence of time dependence of mean 
square displacement (up to 5 x 10 time steps) . The time 
dependencies of temperature, pressure and configuration 
energy were additionally analyzed to control equilibra- 
tion. Data were subsequently collected during the time 
^samp that was chosen to be large enough for correct cal- 
culation of diffusion coefficients by Einstein relation. So 
this time was approximately equal to tgamp ^ 3r, where 
T is the time necessary to reach diffusive regime after 
ballistic and plateau ones. For more details see Supple- 
mentary Material. 

Avoiding spontaneous crystallization during equilibra- 
tion process is the principal difficulty of MD simulations 
of glassy dynamics of liquids. For model glass-forming 
systems this problem is usually solved by either using 
of non-isotropic potentials [5] [52], or considering mul- 
ticomponent systems [23l [24], or using nonequilibrium 
cooling [iHl [251 . For collapsing soft spheres system, it 
is possible to avoid crystallization in the one-component 
system with isotropic potential due to the quasibinary 
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FIG. 3. (Color online) (a) Tg(p) [using (VF)]. (b) Quasibi- 
nary index, F, for T = 0.1 and the range of densities where 
the glassy dynamics was detected. Comparing the F{p) with 
Tg{p) we see that these dependencies are opposite to one an- 
other: the bigger F the smaller Tg and the maximum of the 
former curve is located at the same density as the minimum 
of the later one. The inset shows the radial distribution func- 
tion g{r). It has two clear peaks a.t r — d, r = a that proves 
the quasibinary behaviour of the system, (c) The diffusion 
coefficient D{p) for T = 0.1. 



behavior that develops itself in certain density interval. 
In our case this range is p € (0.51 — 0.74). In the inset in 
Fig.jSj^ we show the splitting of the first peak in the radial 
distribution functions (RDF) that illustrates the quasi- 
binary behaviour of our system in hand. Outside this 
interval, it is hardly possible to study the glassy dynam- 
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FIG. 4. (Color online) (a) Mean-square displacement, and (b) 
Intermediate scattering function for the density p = 0.6. In- 
set shows the snapshot of the typical distribution of particles 
for T = 0.1 and p = 0.6. We mark by red (blue) color parti- 
cles that have a neighbor at distance of the order of the first 
(second) maximum in rdf [hard-core and soft-core diameters] . 



ics because system spontaneously crystalizes whereupon 
supercooling below the melting line. Conversely, inside 
the region mentioned, one can equilibrate supercooled 
liquid without crystallization down to temperatures at 
which relaxation time becomes too large for simulation. 
In our case, these minimal temperatures were chosen so 
that diffusion coefficients were on the order of 10~^. At 
that temperatures, the total calculation time required for 
equilibration and data collecting is ~ 3 x 10^ MD steps 
[three days of calculation on 32 processors in parallel] . 

In order to control stability of glassy state, we per- 
formed calculation at lower temperatures also (down to 
Tg). At that temperatures, the system cannot be equili- 
brated completely because of large relaxation times and 
so we did not use these data for Tg calculations. But 
what we have observed is the absence of any crystalliza- 
tion up to 10^ MD steps and so we conclude that the 
glass state is stable (at least at simulation time scales). 

To acquire information about the glass state we fo- 
cus on the temperature range T > Tg where the "glass- 
forming fluctuations" |27| slow down the system dynamics 
with temperature decreasing. The conventional correla- 



tion function tools have been used: the mean square dis- 
placement Ar^(i) and the intermediate scattering func- 
tion Fs{q,t). The time dependencies of these func- 
tions for p — 0.6 and different temperatures are shown 
in Fig. [4] One can see the typical picture glassform- 
ers demonstrate at low temperature [28j. Namely, the 
"plateau" reflecting the cage effect (when the particle is 
trapped in the "cage" of the nearest neighbours) appears 
on Ar^(t) and Fs{q,t) at sufhciently low temperatures 
that indicates the onset of glassy regime in system dy- 
namics. Meanwhile, the system remains in disordered 
state as can be seen from radial distribution functions, 
see inset of Fig. [Sfj. We note the splitting of the first and 
second peaks of the radial distribution function (Fig. |3]d) . 
The splitting of the first peak reflects the quasibinarity of 
the system caused by the form of the potential (see dis- 
cussion below) . While the splitting of the second peak is 
apparently the (system independent) attribute of glassy 
state [3 [HI US]- 

In order to estimate the glass transition temperature, 
we calculated diffusion coefficient D at each of the tem- 
perature and density investigated. The temperature de- 
pendencies D{T) were approximated by both the Vogel- 
Fulcher formula (VF) [Ml [30] 



D 



— ^0 



exp 



A 

T-T 



(vf) 



(2) 



and the power law from mode-coupling theory 

(MCT) [301 [SB 



D = D, 



(mc) 



The parameters D^J^\ A, T^^^'' 



(mc) rTi(mc) 



(3) 
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obtained using method of nonlinear least squares. Since 
the expressions ([2|, (|3| are correct only in a vicinity of 
glass transition temperature, the temperature interval 
T e (Tlnin, ^max) for least squares approximation was 
chosen so that D{T^^^) ~ 10-^ D{T^^^) ~ lO'^ that 
approximately corresponds T/Tg G (1.15, 1.6). The typi- 
cal temperature dependence D{T) of diffusion coefficient 
obtained from simulation and its VF and MCT approxi- 
mations are shown in the inset of Fig. [2| One can see that 
both formulas provide good fitting of simulation data. 

Having the parameters of ([2]), ([3| one can get the glass 
transition temperature. According to generally accepted 
definition, glass transition occurs at Dq/D = 10", where 
13 < n- < 17. Using ([2]), (|3]) we obtain 



rp{vf) ^ y(vf) 



A log 



10' 



(4) 
(5) 



The graph ofTg{p) is shown in Fig.[2]on the phase dia- 
gram (previously obtained in [12]). It follows that Tg{p) 
dependencies obtained using VF and MCT approaches 
are in good agreement with each other. Particularly, the 
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FIG. 5. (Color online) The pair potential U{r) and the split- 
ted first peak of the radial distribution functions. 



both curves are located under the melting line and have 
a minima in the vicinity of density p = 0.6. 

We calculate the fragility index m showing the de- 
viation of D(T) from Arrhenius law and allowing to 
identify the type of glass-forming system according to 
strong-fragile classification. Starting with the definition 
of fragility in the form [321 - 135] 



m 



91ogio (Do/D) 



d{TjT) 
and using Eqs. ([2]),(|4| we get 



(6) 



T=T„ 



(vf) 



A 



-nlnlO + 1 



(7) 



The equation Q particularly shows that the fragility in- 
dex has minimal value rrimin = n corresponding to the 
limit T^"^'^ /A ->■ that turns ^ to Arrhenius law. Fig.[l] 
shows the density dependence of the reduced fragility, 
77i/TO,nin, of our System in comparison to several glass- 
formers. One can see that the dynamics of our system 
is extremely fragile since the mean value of its fragility 
exceeds the one for decalin — one of the most fragile sys- 
tem [14]. Note that m{p) is nonmonotonic and reveals 
clear minimum at p = 0.6 as well as for Tg{p), see Fig. [5] 
Recently it was shown that the increase of the interaction 
softness can lead to the increase of the fragility [M] [3S] . 
The "softness" of our pair potential is large at r corre- 
sponding to the peaks of the radial distribution function, 
see Fig. [5] Away from p* = 0.6 one of the RDF peaks 
dominates, see inset in Fig. [3]d, that effectively makes 
softer the effective interparticle interaction and helps to 
interpret the nonmonotonic behaviour of m[p). 

Advances over the last decade have linked non- 
Arrhenius behavior of fragile glass formers with the pres- 
ence of locally heterogeneous dynamics: i.e. the presence 
of distinct (if transient) slow and fast regions within the 
material ^136j . The quasibinary character of our particle 



system allows to make the selection among the particles: 
we mark by red (blue) color particles that have a neighbor 
at distance of the order of the first (second) maximum 
in rdf [hard-core and soft-core diameters] , see Figjsj The 
snapshot of the spatial particle distribution over the sim- 
ulation volume (see inset in Fig. |4] and Fig. 2 in the Sup- 
plementary) shows the high degree of the heterogeneity in 
our system that favors the existence of the locally hetero- 
geneous dynamics within the system because of different 
free volumes (and so diffusion coefficients) of particles 
with different effective radii. 

It has been mentioned above that the system demon- 
strates quasibinary behavior due to the repulsive shoul- 
der of the pair potential, see inset in Fig. [3]d. It refiects 
the competition between hard-core and soft-core scales 
and, as a result, the frustration in the system. As the 
quasibinarity index we choose 



F = 



2piP2 

pI + pV 



(8) 



where pi 



p f^^ r^a(r)dr P I^^ r^g(r)dr 

};i r- J ^i^d P2 = r2dr " Here n 

and r2 are minima of the first and the second RDF peaks. 
So pi , p2 are local densities in the vicinity of the peaks 
and F is their inverse symmetrized ratio. 

The origin of the frustration in our system is the same 
as for binary soft-sphere like systems where the frustra- 
tion is due to the large, local rearrangement of atoms re- 
quired for the formation of a crystal from a fluid or glassy 
configuration f37'. This situation is caused by presence 
of second component and so the frustration of such type 
increases with increasing concentration. Thus the quasi- 
binary index F can serve as frustration measure in our 
system. It is clear that the value ~ corresponds to 
the situation than only one of the scales dominates and 
so there is no frustration in the system. On the contrary, 
if F ~ 1 the competition between hard-core and soft- 
core scales and so the frustration are maximal. It should 
be noted that F depends on temperature as well as on 
density. In order to study the influence of frustration on 
glass-forming ability of the system we calculate the den- 
sity dependence of the frustration index at sufficiently 
low temperatures. In Fig. [sja we show F{p) at T = 0.08 
(of the order of Tn for all the densities investigated). It 
follows from Fig.lSlthat F{p) and Tg{p) have opposite be- 
haviour: the bigger F the smaller Tg and the maximum 
of the former curve is located at the same density as the 
minimum of the later one. 

At the density p* where Tg (p) has minimum the system 
has potentially a number of different lattice constants as 
we mentioned above. If we imagined the long range order 
formation at p* then the crystal would be strongly dis- 
torted by defects due to the competition of the different 
lattice constants. This situation should favour the diffu- 
sion and frustration at the same time, see Fig. |3] Tg is 
determined both by the diffusion and by the frustration: 
diffusion tries to decrease Tg while the frustration does 
the opposite. However the diffusion defeats frustration 
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in our system, so Tg has minimum at p* . 

Finally we discuss the boundaries of the density do- 
main, p £ (0.35 — 0.75). Beyond these boundaries the 
glass forming ability quickly decreases. This one can 
judge from the quick rise of the diffusion coefficient, D{p), 
at the melting line, see inset in Fig. [3]:, and the destruc- 
tion of the quasibinary behaviour, see inset in Fig. [3}d. 

In conclusion, we show using the molecular dynam- 
ics that the one component simple liquid with pure re- 
pulsive potential shows glassy behaviour in quasiequilib- 
rium cooling. We explain the nonmonotonic density de- 
pendence of Tg, frustration, the diffusion coefficient and 
fragility by the evolution of the quasibinary properties 
of our system. We observe that our system belongs to 
the short list of the most fragile systems. Our system 
can be used as the simple toy model for investigation 
of the quasibinary frustrated systems. In the search for 
(super)fragility in the one component repulsive simple 
liquids one should test them for quasibinarity. 
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SUPPLEMENTARY MATERIAL 

Appendix A: Details of simulation and results 
handling 

1. MD simulation 

For MD simulations, we have used DL_POLY 2.20 
Molecular Simulation Package 01] developed at Dares- 
bury Laboratory. 

Initially, the system was equilibrated from lattice con- 
figuration at high temperature for 10^ MD time steps 
and then it was consecutively equilibrated at lower tem- 
peratures from final configurations of previous runs. At 
relatively high temperatures (T > 1.2rg), the equilibra- 
tion time was chosen to be in order of 10'* — 2 x 10^ time 
steps. The time dependencies of temperature, pressure 
and configuration energy were analyzed to control equili- 
bration. At low (T < 1.2rg) temperatures addition con- 
trol is needed so systems was consecutively equilibrated 
until convergence of time dependence of mean square dis- 
placement. Data were subsequently collected during the 
time tsamp that was chosen to be large enough for cor- 
rect calculation of diffusion coefficients by Einstein rela- 
tion. So this time was approximately equal to tsamp ^ 3t, 
where r is the relaxation time, see Fig.[6]^a). The reason 
is simple: one r is needed to reach diffusive regime after 



ballistic and plateau ones, another r is the main time in- 
terval for calculation by Einstein formula and one more 
T is to improve averaging statistics for MSD calculations 
in the main interval. 

We obtain Tg by quasiequilibrium method that is 
just approximation of equilibrium D{T) curve for super- 
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FIG. 6. (Color online) (a)"3r rule" for correct calculation of 
diffusion coefficients from mean square displacement. This 
rule says that the total time for data collecting should be in 
order of 3r, where t is the relaxation time. The reason is 
simple: one r is needed to reach diffusive regime after ballis- 
tic and plateau ones, another r is the main time interval for 
calculation by Einstein formula and one more r is to improve 
averaging statistics for MSD calculations in the main inter- 
val, (b) The potential of the collapsing soft spheres model, (c) 
Angel plot for temperature dependence of diffusion coefficient 
demonstrating essential deviation of D(T) from Arrhenius law 
and so the fragile character of system dynamics. 
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cooled liquid by appropriate fitting formula (for example, 
Vogel-Fulcher law) and using the relation D{Tg)/Do = 
10", where n = 13 — 17. This method is natural for 
(quasi)equilibrium simulation we have used but it differs 
from the dynamic method. According to the dynamic 
method, glass transition temperature is obtained as tem- 
perature at which the system falls out of equilibrium (i.e. 
when the cooling rate exceeds the relaxation one). This 
temperature can be determined by sharp peaks on heat 
capacity Cp(T) or thermal expansion coefficient, ap{T). 
The temperature, Tg, obtained by the dynamic method, 
depends on the coofing rate v, i.e. Tg — Tg{v). The 
natural question appears: how do the glass transition 
temperatures obtained by static (quasiequilibrium) and 
dynamic methods correlate with each other? The an- 
swer is not conclusive because the rigorous theory is still 
absent but some arguments may be given. There are a 
number of experimental and theoretical works (see, e.g., 
Refs. [2, 4, 25, 39-43^) where Tg{v) dependence has been 
investigated. It was shown that this dependence can be 
represented in the form. 



Tg{v)=Tg{V^O)+fiv), 



(Al) 



where f{v) is some function of the cooling rate v = 
where T and t are temperature and time correspondingly. 
The function f{v) is often approximated as follows: 



fiv) 



B 

In {vo/v) ■ 



(A2) 



where B and vq are constants [39l |4T] . Expression ( A2 ) 



is able to give a satisfactory fit to the experimental data 
when V is varied over 3 decades (Ml. The shape of 



Eq. ( A2 ) can be justified if we use the Vogel-Fulcher ap- 



proximation for the relaxation time. 



t(T) = Tocxp 



B 



T~Tg{v^Q) 



(A3) 



and if we assume [see, e.g., Refs. [321 SD US]] that the 
system falls out of equilibrium at temperature Tg{v) 
where the relaxation time is of the order of the inverse 
of the cooling rate, i.e. Tg{v)/T\Tg{vy\ ^ v [and so 



Eq. (A2) is applicable while v < 



It follows from the derivation that 
Wo, see, e.g., Ref. [33]. 
Note that there are other suggestions for f(v) that also 
agree not so bad with experiments 135] . In any 

case, at sufficiently small cooling rate Tg becomes nearly 
independent on the cooling rate v. This is the limit where 
we perform our simulations. Then pressure and temper- 
ature are well defined in Nozier-Hover thermostat that 
we have used. 



corresponding to hard- and soft-cores, see Fig.[6];b). This 
potential was introduced in Refs. [12l (TS] [16] in the fol- 
lowing form: 

U{r)^s(^^y +'-{l-ts.nh[ko{r-a,)]}. (A4) 

Using the identity 

nF(a;) = ^{l-tanh|}, (A5) 

where the Fermi- function, np^x) = 1/[1 -|- exp(x)], we 
rewrote the potential of the collapsing soft spheres model: 

U{r)=ei^-j +enp[2ko{r-<j^)]. (A6) 

3. Angell plot 

Angell plot for temperature dependence of diffusion co- 
efficient demonstrating essential deviation of D(T) from 
Arrhenius law and so the fragile character of system dy- 
namics, see Fig. [gJc). The construction of such plot is 
the traditional way to distinguish between the strong and 
fragile glassformer [13]. 

Appendix B: Discussion 

What is the basis for the paradigm that no glass can 
be expected for purely repulsive potential for quasiequi- 
librium cooling? First of all this is "experimental" fact. 




2. The potential of the collapsing soft spheres 
model 

The potential of the collapsing soft spheres model is 
the pure repulsive potential with two characteristic scales 



FIG. 7. (Color online) Snapshot of the typical distribution 
of particles over the simulation volume. Here T = 0.1, and 
p = 0.6. We mark by red (blue) color particles that have a 
neighbor at distance « ri (r2) [ri and r2 are defined in Fig. 5 
of our paper] . 
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Namely, it is definitely established last decades from nu- 
merical molecular-dynamic experiments that purely re- 
pulsive potential in one-component isotropic system does 
not allow to freeze the system into a glass while cooling 
is quasiequilibrium (as it was defined above): the sys- 
tem spontaneously crystalizes (there is lack of frustration 
built in this system), see also discussion in Refs. 
However, if the cooling rate is large enough, larger than 
the typical relaxation rate, then it is possible to freeze 
such system into a glass with Tg essentially depending on 
the cooling rate. But in this dynamical cooling regime 
the region above Tg, we focus on in our paper, is not 
accessible where glassy slowing down is already present 
but the system is still in the liquid state. The dynam- 
ical cooling is in fact the method to "cook" the glass 
and find Tg{v). Both regimes, qusiequilibrium and dy- 
namical, are accessible experimentally. Therefore, if we 
want to study the mechanism of glassy slow down in su- 
percooled state we have to use quasiequilibrium method. 
So our model system should not spontaneously crystal- 
lize being supercooled in quasiequilibrium conditions due 
to some frustrations. If we will speak on the intuitive 
level it becomes difficult to imagine frustration in simple 
monatomic pure repulsive system: potential is isotropic 
and what is most important only one characteristic scale 



is built in. So we propose the repulsive potential with 
additional source of frustration. Our potential has the 
soft step and the simple liquid with this interaction has 
two characteristic scales. The competition between these 
scales makes the system effectively quasibinary. And this 
is the origin of frustration in our system. 

The physical origin of the non-Arrhenius behavior of 
fragile glass formers is an area of active investigation in 
glass physics. Advances over the last decade have linked 
this phenomenon with the presence of locally heteroge- 
neous dynamics in fragile glass formers; i.e. the presence 
of distinct (if transient) slow and fast regions within the 
material [2 . The quasibinary character of our particle 
system allows to make the selection among the particles: 
we mark by red (blue) color particles that have a neighbor 
at distance of the order of the first (second) maximum 
in rdf [hard-core and soft-core diameters], see Fig. 5 in 
the main paper. To demonstrate the high degree of the 
heterogeneity in our system we place the snapshot of the 
spatial particle distribution over the simulation volume, 
see Fig. [7] The high spatial heterogeneity of our system 
favors the existence of the locally heterogeneous dynam- 
ics within the system because of different free volumes 
(and so diffusion coefficients) of particles with different 
effective radii. 
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